execute mcmc sampling

goMCMC=function(mdl,data,smp=500,wrm=100,th=1){
  mcmc=mdl$sample(
  data=data,
  seed=1,
  chains=4,
  iter_sampling=smp,
  iter_warmup=wrm,
  thin=th,
  refresh=1000
  )
  mcmc
}

see mcmc result and parameters

seeMCMC=function(mcmc,exc='',ch=T){ # exclude 'exc' parameters from seeing
  print(mcmc)
  prs=mcmc$metadata()$model_params[-1] # reject lp__
  for(pr in prs){
    if(grepl('^y',pr)) next # not show predictive value "y*" information
    if(exc!='' && grepl(paste0('^',exc),pr)) next
    drw=mcmc$draws(pr)
    if(ch){
      par(mfrow=c(1,4))
      drw[,1,] |> plot(type='l',main='chain1',ylab=pr)
      drw[,2,] |> plot(type='l',main='chain2',ylab=pr)
      drw[,3,] |> plot(type='l',main='chain3',ylab=pr)
      drw[,4,] |> plot(type='l',main='chain4',ylab=pr)
    }

    par(mfrow=c(1,2))
    drw |> hist(main=pr,xlab='')
    drw |> density() |> plot(main=pr)    
  }
}

ex1 normal distribution
ex1.stan

data{
  int N;
  vector[N] y;
}

parameters{
  real m;
  real<lower=0> s;
}

model{
  y~normal(m,s);
}
mdl=cmdstan_model('./ex1.stan')
y=rnorm(10,2,1)
summary(y)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.662   1.042   2.136   1.901   2.546   2.904
par(mfrow=c(1,2))
hist(y,main='y')
density(y) |> plot(main='y')  

data=list(N=length(y),y=y)

mle=mdl$optimize(data=data)
## Initial log joint probability = -1357.56 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       12       -3.2695   4.48716e-05   1.26718e-06           1           1       23    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.3 seconds.
mle
##  variable estimate
##      lp__    -3.27
##      m        1.90
##      s        0.84
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.3 seconds.
mcmc
##  variable  mean median   sd  mad    q5   q95 rhat ess_bulk ess_tail
##      lp__ -4.51  -4.17 1.15 0.78 -6.75 -3.46 1.00      802      856
##      m     1.93   1.93 0.34 0.30  1.39  2.48 1.01      890      693
##      s     1.04   0.98 0.32 0.24  0.68  1.60 1.00     1194     1099
mcmc$metadata()$stan_variables
## [1] "lp__" "m"    "s"
mcmc$metadata()$model_params
## [1] "lp__" "m"    "s"
seeMCMC(mcmc)
##  variable  mean median   sd  mad    q5   q95 rhat ess_bulk ess_tail
##      lp__ -4.51  -4.17 1.15 0.78 -6.75 -3.46 1.00      802      856
##      m     1.93   1.93 0.34 0.30  1.39  2.48 1.01      890      693
##      s     1.04   0.98 0.32 0.24  0.68  1.60 1.00     1194     1099

ex2 poisson distribution

ex2

y=rpois(20,3)
summary(y)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00    1.00    2.50    2.65    4.00    6.00
par(mfrow=c(1,2))
hist(y,main='y')
density(y) |> plot(main='y')

data=list(N=length(y),y=y)

ex2-1.stan

data{
  int N;
  vector[N] y;
}

parameters{
  real<lower=0> m;
  real<lower=0> s;
}

model{
  y~normal(m,s);
}

generated quantities{
  vector[N] y1;
  for(i in 1:N)
    y1[i]=normal_rng(m,s);
}
# when fitting poisson dist. sample to normal dist.
mdl=cmdstan_model('./ex2-1.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -38.2607 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##        7      -20.3939   0.000406511   0.000836652           1           1        8    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##    lp__     -20.39
##    m          2.65
##    s          1.68
##    y1[1]      3.39
##    y1[2]      0.92
##    y1[3]     -0.01
##    y1[4]      3.33
##    y1[5]      1.89
##    y1[6]      2.91
##    y1[7]      3.97
##    y1[8]      2.69
##    y1[9]      2.26
##    y1[10]     2.76
##    y1[11]     4.05
##    y1[12]     3.32
##    y1[13]     5.38
##    y1[14]     3.00
##    y1[15]     3.31
##    y1[16]     6.28
##    y1[17]     2.42
##    y1[18]     4.76
##    y1[19]     1.93
##    y1[20]    -1.42
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -19.95 -19.60 1.04 0.76 -22.10 -18.94 1.00      757     1110
##    m        2.66   2.66 0.41 0.41   1.99   3.33 1.01     1420     1128
##    s        1.85   1.80 0.33 0.31   1.39   2.49 1.00     1180      952
##    y1[1]    2.66   2.64 1.95 1.81  -0.58   5.83 1.00     1658     1904
##    y1[2]    2.63   2.70 1.97 1.84  -0.62   5.85 1.00     1771     1317
##    y1[3]    2.66   2.65 1.94 1.90  -0.49   5.84 1.00     2050     1747
##    y1[4]    2.67   2.63 1.88 1.82  -0.36   5.72 1.00     2008     2059
##    y1[5]    2.62   2.59 1.93 1.88  -0.48   5.82 1.00     1916     1822
##    y1[6]    2.68   2.66 1.94 1.85  -0.50   5.92 1.00     1980     1999
##    y1[7]    2.66   2.71 1.88 1.76  -0.40   5.82 1.00     2004     1933
##    y1[8]    2.69   2.61 1.90 1.80  -0.28   5.94 1.00     1867     2007
##    y1[9]    2.63   2.63 1.94 1.83  -0.54   5.83 1.00     1984     1969
##    y1[10]   2.70   2.73 1.89 1.74  -0.51   5.77 1.00     1995     1876
##    y1[11]   2.71   2.70 1.87 1.84  -0.35   5.82 1.00     2019     1854
##    y1[12]   2.62   2.67 1.86 1.87  -0.44   5.61 1.00     2097     1854
##    y1[13]   2.60   2.60 1.92 1.82  -0.53   5.86 1.00     1924     1901
##    y1[14]   2.68   2.66 1.89 1.81  -0.42   5.73 1.00     2003     1938
##    y1[15]   2.71   2.75 1.98 1.96  -0.53   6.10 1.00     2085     1972
##    y1[16]   2.59   2.62 1.95 1.89  -0.56   5.79 1.00     1896     1974
##    y1[17]   2.65   2.58 1.91 1.85  -0.51   5.83 1.00     1958     2052
##    y1[18]   2.66   2.68 1.94 1.90  -0.47   5.71 1.00     2017     1919
##    y1[19]   2.65   2.60 1.94 1.87  -0.58   5.89 1.00     1902     1954
##    y1[20]   2.66   2.68 1.97 1.84  -0.69   5.91 1.00     1950     1971

y1=mcmc$draws('y1')
par(mfrow=c(1,2))
hist(y1,main='y1')
density(y1) |> plot(main='y1')     

ex2-2.stan

data{
  int N;
  array[N] int y;
}

parameters{
  real<lower=0> l;
}

model{
  y~poisson(l);
}

generated quantities{
  array[N] int y1;
  for(i in 1:N)
    y1[i]=poisson_rng(l);
}
mdl=cmdstan_model('./ex2-2.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -1.52152 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##        4      -1.34834   1.36669e-05   7.81676e-06        0.99        0.99        6    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##    lp__      -1.35
##    l          2.65
##    y1[1]      4.00
##    y1[2]      2.00
##    y1[3]      3.00
##    y1[4]      2.00
##    y1[5]      2.00
##    y1[6]      0.00
##    y1[7]      2.00
##    y1[8]      2.00
##    y1[9]      3.00
##    y1[10]     1.00
##    y1[11]     2.00
##    y1[12]     2.00
##    y1[13]     5.00
##    y1[14]     2.00
##    y1[15]     0.00
##    y1[16]     2.00
##    y1[17]     5.00
##    y1[18]     5.00
##    y1[19]     1.00
##    y1[20]     3.00
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc)
##  variable  mean median   sd  mad    q5   q95 rhat ess_bulk ess_tail
##    lp__   -0.84  -0.57 0.68 0.28 -2.16 -0.37 1.00     1007     1383
##    l       2.70   2.68 0.36 0.35  2.14  3.32 1.00      812     1179
##    y1[1]   2.70   3.00 1.66 1.48  0.00  6.00 1.00     1813     1763
##    y1[2]   2.67   2.00 1.69 1.48  0.00  6.00 1.00     1971     1892
##    y1[3]   2.73   3.00 1.69 1.48  0.00  6.00 1.00     1894     1902
##    y1[4]   2.68   2.00 1.64 1.48  0.00  6.00 1.00     1890     1918
##    y1[5]   2.67   2.00 1.69 1.48  0.00  6.00 1.00     1933     2122
##    y1[6]   2.72   3.00 1.69 1.48  0.00  6.00 1.00     1744     1890
##    y1[7]   2.66   3.00 1.61 1.48  0.00  5.00 1.00     2048     2033
##    y1[8]   2.70   3.00 1.70 1.48  0.00  6.00 1.00     2075     1933
##    y1[9]   2.73   3.00 1.71 1.48  0.00  6.00 1.00     1883     1978
##    y1[10]  2.67   2.00 1.65 1.48  0.00  6.00 1.00     1697     1875
##    y1[11]  2.67   3.00 1.65 1.48  0.00  6.00 1.00     1952     1916
##    y1[12]  2.66   2.00 1.67 1.48  0.00  6.00 1.00     1851     1833
##    y1[13]  2.68   3.00 1.70 1.48  0.00  6.00 1.00     1886     1936
##    y1[14]  2.74   3.00 1.66 1.48  0.00  6.00 1.00     1888     1829
##    y1[15]  2.70   2.00 1.73 1.48  0.00  6.00 1.00     1923     1924
##    y1[16]  2.64   2.00 1.66 1.48  0.00  6.00 1.00     1797     1824
##    y1[17]  2.75   3.00 1.68 1.48  0.00  6.00 1.00     1710     1770
##    y1[18]  2.66   2.00 1.66 1.48  0.00  6.00 1.00     2051     1845
##    y1[19]  2.70   2.00 1.70 1.48  0.00  6.00 1.00     1799     2001
##    y1[20]  2.65   2.00 1.64 1.48  0.00  6.00 1.00     1986     2102

y1=mcmc$draws('y1')
par(mfrow=c(1,2))
hist(y1,main='y1')
density(y1) |> plot(main='y1') 

ex3
difference test

ex3.stan

data{
  int N;
  vector[N] a;
  vector[N] b;
}
parameters{
  real ma;
  real<lower=0> sa;
  real mb;
  real<lower=0> sb;
}
model{
  a~normal(ma,sa);
  b~normal(mb,sb);
}
generated quantities{
  real d;
  d=mb-ma;
}
n=20
a=rnorm(n,10,1)
b=rnorm(n,11,2)
data=list(N=n,a=a,b=b)

mdl=cmdstan_model('./ex3.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -6122.02 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       29      -35.2022   0.000122699   0.000172573           1           1       44    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##      lp__   -35.20
##      ma      10.16
##      sa       1.10
##      mb      11.21
##      sb       1.94
##      d        1.05
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,ch=F)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##      lp__ -36.45 -36.15 1.42 1.27 -39.16 -34.79 1.00      792     1157
##      ma    10.17  10.16 0.27 0.25   9.73  10.61 1.00     1388     1303
##      sa     1.21   1.18 0.21 0.20   0.92   1.61 1.00     2003     1484
##      mb    11.21  11.22 0.47 0.46  10.41  11.96 1.01      822      833
##      sb     2.14   2.09 0.38 0.35   1.61   2.84 1.00     2235     1408
##      d      1.04   1.04 0.54 0.53   0.14   1.92 1.00      955      887

d=mcmc$draws('d')
d
## # A draws_array: 500 iterations, 4 chains, and 1 variables
## , , variable = d
## 
##          chain
## iteration    1    2    3    4
##         1 0.94 0.50 0.96 1.51
##         2 1.27 0.69 0.77 1.43
##         3 1.44 0.63 0.39 0.45
##         4 1.44 0.55 0.64 2.80
##         5 1.54 1.75 0.76 2.04
## 
## # ... with 495 more iterations
mean(d>0)
## [1] 0.971

ex4

ex4-1.stan

data{
  int N;
  vector[N] x;
  vector[N] y;
}

parameters{
  real b0;
  real b1;
  real<lower=0> s;
}

model{
  y~normal(b0+b1*x,s);
}
# make prediction with R
n=20
x=runif(n,0,100)
y=rnorm(n,x*3+10,10)
par(mfrow=c(1,1))
plot(x,y)

data=list(N=n,x=x,y=y)

mdl=cmdstan_model('./ex4-1.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -291608 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       36      -50.5368    0.00179027   0.000612657      0.9903      0.9903       84    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##      lp__   -50.54
##      b0       8.72
##      b1       3.00
##      s        7.59
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.2 seconds.
## Chain 2 finished in 0.2 seconds.
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 finished in 0.2 seconds.
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 finished in 0.2 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.2 seconds.
## Total execution time: 0.3 seconds.
seeMCMC(mcmc,ch=F)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##      lp__ -50.03 -49.72 1.21 1.02 -52.31 -48.68 1.00      510      474
##      b0     8.68   8.58 4.16 3.92   1.96  15.52 1.01      249      260
##      b1     3.00   3.00 0.07 0.06   2.89   3.11 1.01      305      375
##      s      8.63   8.47 1.57 1.49   6.45  11.37 1.00      931     1045

b0=mcmc$draws('b0')
b1=mcmc$draws('b1')
b=tibble(b0=c(b0),b1=c(b1)) |> sample_n(100)
x1=runif(nrow(b),min(x),max(x))
xy=tibble(x=x1,y=b$b0+b$b1*x1)
par(mfrow=c(1,1))
plot(xy)

ex4-2.stan

data{
  int N;
  vector[N] x;
  vector[N] y;
  int N1;
  vector[N1] x1;
}

parameters{
  real b0;
  real b1;
  real<lower=0> s;
}

model{
  y~normal(b0+b1*x,s);
}

generated quantities{
  vector[N1] m;
  vector[N1] y1;
  for(i in 1:N1){
    m[i]=b0+b1*x1[i];
    y1[i]=normal_rng(m[i],s);
  }
}
# make prediction with stan
n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition

data=list(N=n,x=x,y=y,N1=n1,x1=x1)

mdl=cmdstan_model('./ex4-2.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -3431.79 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       30      -50.5368     0.0262534   0.000427242      0.9908      0.9908       54    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##    lp__     -50.54
##    b0         8.72
##    b1         3.00
##    s          7.59
##    m[1]      48.61
##    m[2]      76.03
##    m[3]     103.44
##    m[4]     130.86
##    m[5]     158.27
##    m[6]     185.69
##    m[7]     213.10
##    m[8]     240.52
##    m[9]     267.93
##    m[10]    295.35
##    y1[1]     44.85
##    y1[2]     80.67
##    y1[3]    102.04
##    y1[4]    131.69
##    y1[5]    150.33
##    y1[6]    190.10
##    y1[7]    215.18
##    y1[8]    232.94
##    y1[9]    272.21
##    y1[10]   301.88
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.2 seconds.
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 finished in 0.2 seconds.
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 finished in 0.2 seconds.
## Chain 4 finished in 0.2 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.2 seconds.
## Total execution time: 0.4 seconds.
mcmc$metadata()$stan_variables
## [1] "lp__" "b0"   "b1"   "s"    "m"    "y1"
seeMCMC(mcmc,'m',ch=F)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -50.09 -49.75 1.30 1.03 -52.55 -48.71 1.01      421      603
##    b0       9.23   9.18 4.42 4.31   1.95  16.95 1.01      216      258
##    b1       2.99   2.99 0.07 0.07   2.87   3.10 1.01      256      379
##    s        8.61   8.39 1.50 1.42   6.57  11.28 1.00      885      798
##    m[1]    49.03  48.89 3.61 3.48  43.02  55.37 1.01      219      265
##    m[2]    76.38  76.24 3.09 2.95  71.24  81.70 1.01      227      274
##    m[3]   103.73 103.63 2.62 2.58  99.56 108.13 1.01      248      302
##    m[4]   131.09 131.06 2.25 2.22 127.62 134.66 1.01      303      309
##    m[5]   158.44 158.44 2.01 1.96 155.24 161.59 1.00      476      547
##    m[6]   185.79 185.82 1.96 1.86 182.63 188.85 1.00     1154     1260
##    m[7]   213.15 213.12 2.11 2.01 209.77 216.52 1.00     2077     1397
##    m[8]   240.50 240.54 2.43 2.36 236.65 244.47 1.00     1758     1322
##    m[9]   267.85 267.83 2.86 2.77 263.29 272.65 1.00     1006     1329
##    m[10]  295.21 295.17 3.36 3.23 289.88 300.81 1.00      699     1185
##    y1[1]   49.35  49.14 9.48 9.42  34.30  64.92 1.00     1065     1494
##    y1[2]   76.71  76.80 9.06 8.75  61.65  91.33 1.00     1288     1466
##    y1[3]  103.47 103.65 9.11 8.83  88.77 118.28 1.00     1355     1780
##    y1[4]  130.97 131.04 8.78 8.47 116.64 145.04 1.00     1730     1730
##    y1[5]  159.18 159.02 9.04 8.59 143.77 174.21 1.00     1909     1889
##    y1[6]  185.55 185.52 8.91 8.76 171.16 200.23 1.00     2092     1887
##    y1[7]  213.14 213.06 8.90 8.75 198.91 227.82 1.00     1653     1816
##    y1[8]  240.40 240.56 9.16 8.88 225.17 255.09 1.00     1820     1894
##    y1[9]  268.00 267.95 9.33 9.00 252.95 283.27 1.00     1960     1932
##    y1[10] 295.37 295.19 9.37 8.58 280.13 310.57 1.00     1565     1649

m=mcmc$draws('m')
summary(m)
## # A tibble: 10 × 10
##    variable  mean median    sd   mad    q5   q95  rhat ess_bulk ess_tail
##    <chr>    <dbl>  <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>    <dbl>    <dbl>
##  1 m[1]      49.0   48.9  3.61  3.48  43.0  55.4  1.01     220.     265.
##  2 m[2]      76.4   76.2  3.09  2.95  71.2  81.7  1.01     228.     274.
##  3 m[3]     104.   104.   2.62  2.58  99.6 108.   1.01     248.     302.
##  4 m[4]     131.   131.   2.25  2.22 128.  135.   1.01     304.     310.
##  5 m[5]     158.   158.   2.01  1.96 155.  162.   1.00     476.     548.
##  6 m[6]     186.   186.   1.96  1.86 183.  189.   1.00    1155.    1260.
##  7 m[7]     213.   213.   2.11  2.01 210.  217.   1.00    2078.    1398.
##  8 m[8]     240.   241.   2.43  2.36 237.  244.   1.00    1758.    1322.
##  9 m[9]     268.   268.   2.86  2.77 263.  273.   1.00    1006.    1330.
## 10 m[10]    295.   295.   3.36  3.23 290.  301.   1.00     700.    1186.
smm=summary(m)

y1=mcmc$draws('y1')
summary(y1)
## # A tibble: 10 × 10
##    variable  mean median    sd   mad    q5   q95  rhat ess_bulk ess_tail
##    <chr>    <dbl>  <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>    <dbl>    <dbl>
##  1 y1[1]     49.3   49.1  9.48  9.42  34.3  64.9  1.00    1065.    1494.
##  2 y1[2]     76.7   76.8  9.06  8.75  61.7  91.3  1.00    1289.    1466.
##  3 y1[3]    103.   104.   9.11  8.83  88.8 118.   1.00    1356.    1780.
##  4 y1[4]    131.   131.   8.78  8.47 117.  145.   1.00    1731.    1730.
##  5 y1[5]    159.   159.   9.04  8.59 144.  174.   1.00    1909.    1890.
##  6 y1[6]    186.   186.   8.91  8.76 171.  200.   1.00    2093.    1888.
##  7 y1[7]    213.   213.   8.90  8.75 199.  228.   1.00    1653.    1817.
##  8 y1[8]    240.   241.   9.16  8.88 225.  255.   1.00    1821.    1894.
##  9 y1[9]    268.   268.   9.33  9.00 253.  283.   1.00    1961.    1932.
## 10 y1[10]   295.   295.   9.37  8.58 280.  311.   1.00    1565.    1650.
smy=summary(y1)

xy=tibble(x=x1,m=smm$median,ml5=smm$q5,mu5=smm$q95,yl5=smy$q5,yu5=smy$q95)

par(mfrow=c(1,1))
xlim=c(min(x),max(x))
ylim=c(min(y),max(y))
plot(x,y,
     xlim=xlim,ylim=ylim, xlab='x',ylab='y',col='red')
par(new=T)
plot(xy$x,xy$m,type='l',
     xlim=xlim,ylim=ylim, xlab='',ylab='',col='black')
par(new=T)
plot(xy$x,xy$ml5,type='l',
     xlim=xlim,ylim=ylim, xlab='',ylab='',col='darkgray')
par(new=T)
plot(xy$x,xy$mu5,type='l',
     xlim=xlim,ylim=ylim, xlab='',ylab='',col='darkgray')
par(new=T)
plot(xy$x,xy$yl5,type='l',
     xlim=xlim,ylim=ylim, xlab='',ylab='',col='lightgray')
par(new=T)
plot(xy$x,xy$yu5,type='l',
     xlim=xlim,ylim=ylim, xlab='',ylab='',col='lightgray')

qplot(x,y,col=I('red'))+
  geom_line(aes(x=x,y=m),xy,col='black')+
  geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
  geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
  geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
  geom_line(aes(x=x,y=yu5),xy,col='lightgray')

ex4-3.stan

data{
  int N;
  vector[N] x;
  vector[N] y;
}

parameters{
  real b0;
  real b1;
  real<lower=0> s;
}

model{
  y~normal(b0+b1*x,s);
}

generated quantities{
  vector[N] m;
  vector[N] y1;
  for(i in 1:N){
    m[i]=b0+b1*x[i];
    y1[i]=normal_rng(m[i],s);
  }
}
data=list(N=n,x=x,y=y)

mdl=cmdstan_model('./ex4-3.stan')

mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.2 seconds.
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 finished in 0.2 seconds.
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 finished in 0.2 seconds.
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 finished in 0.2 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.2 seconds.
## Total execution time: 0.3 seconds.
seeMCMC(mcmc,'m',ch=F)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -50.05 -49.65 1.36 0.99 -52.86 -48.66 1.02      346      218
##    b0       9.02   8.79 4.34 3.79   2.47  16.12 1.05      108      105
##    b1       2.99   2.99 0.07 0.06   2.88   3.10 1.04      153      147
##    s        8.66   8.43 1.62 1.40   6.42  11.64 1.01      667      659
##    m[1]    69.42  69.25 3.14 2.85  64.69  74.50 1.05      102       88
##    m[2]    64.71  64.54 3.23 2.92  59.88  69.98 1.05      103       85
##    m[3]   131.00 130.86 2.16 2.01 127.67 134.50 1.03      150      184
##    m[4]    48.85  48.68 3.54 3.19  43.53  54.59 1.05      106      106
##    m[5]   245.84 245.81 2.38 2.28 241.89 249.83 1.00     1551     1377
##    m[6]   257.86 257.83 2.57 2.50 253.58 262.07 1.00     1186     1137
##    m[7]   293.54 293.53 3.19 3.13 288.23 298.56 1.01      669      826
##    m[8]   144.56 144.42 2.02 1.85 141.43 147.80 1.03      197      260
##    m[9]   256.76 256.73 2.55 2.48 252.50 260.94 1.00     1214     1159
##    m[10]  246.59 246.57 2.40 2.30 242.62 250.59 1.00     1526     1298
##    m[11]  169.93 169.89 1.86 1.71 166.96 173.04 1.01      472      482
##    m[12]  295.20 295.18 3.22 3.16 289.83 300.26 1.01      655      826
##    m[13]  292.87 292.86 3.18 3.12 287.59 297.89 1.01      675      801
##    m[14]  192.33 192.35 1.87 1.72 189.25 195.43 1.00     1104     1300
##    m[15]  112.33 112.19 2.42 2.25 108.73 116.20 1.04      111      135
##    m[16]  272.50 272.47 2.81 2.74 267.74 277.05 1.01      894      978
##    m[17]  128.08 127.95 2.20 2.08 124.71 131.69 1.03      144      178
##    m[18]   94.95  94.80 2.69 2.50  90.94  99.19 1.04      101       91
##    m[19]   91.03  90.90 2.76 2.54  86.90  95.34 1.04      100       91
##    m[20]  115.46 115.33 2.37 2.21 111.90 119.22 1.04      115      135
##    y1[1]   69.53  69.61 9.18 8.84  54.45  84.47 1.01     1087     1936
##    y1[2]   65.05  65.05 9.19 9.05  50.90  80.39 1.01     1024      974
##    y1[3]  131.00 130.90 8.85 8.67 116.54 145.07 1.00     1691     1480
##    y1[4]   48.90  48.73 9.48 8.70  33.72  64.47 1.01      568     1451
##    y1[5]  245.62 245.49 9.17 8.83 230.45 260.66 1.00     2007     1564
##    y1[6]  258.37 258.37 9.11 8.48 243.19 273.56 1.00     1879     2053
##    y1[7]  293.45 293.51 9.40 9.01 278.12 308.61 1.00     1760     1566
##    y1[8]  144.66 144.76 9.00 8.86 129.81 159.51 1.00     1630     1916
##    y1[9]  256.84 257.19 9.18 8.79 240.99 271.41 1.00     1988     1916
##    y1[10] 246.74 246.88 9.13 8.69 231.05 261.44 1.00     1897     1933
##    y1[11] 169.82 169.56 9.22 8.55 154.84 185.25 1.00     2001     1618
##    y1[12] 295.01 295.26 9.21 8.71 279.94 310.23 1.00     1875     1841
##    y1[13] 293.37 293.38 9.13 8.66 278.28 307.77 1.00     1610     1916
##    y1[14] 192.45 192.46 8.90 8.68 178.07 206.77 1.00     2046     1810
##    y1[15] 112.51 112.53 9.14 8.72  98.03 127.82 1.00     1517     1890
##    y1[16] 272.52 272.34 9.28 8.46 257.61 287.76 1.00     1646     1816
##    y1[17] 128.04 128.01 9.12 8.68 113.11 142.53 1.00     1619     1599
##    y1[18]  95.08  95.28 9.06 8.86  80.06 109.78 1.00     1412     1673
##    y1[19]  90.73  90.72 9.20 9.14  75.93 106.04 1.00     1149     1785
##    y1[20] 115.27 115.36 9.18 8.63  99.96 129.84 1.00     1607     1736

y1=mcmc$draws('y1')
smy=summary(y1)

par(mfrow=c(1,1))
plot(y,smy$median,xlab='obs.',ylab='prd.')
abline(0,1)

qplot(y,smy$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1)

par(mfrow=c(1,2))
hist(y-smy$median,xlab='obs.-prd.',main='residual')
density(y-smy$median) |> plot(xlab='obs.-prd.',main='residual')

estimate correlation ex4-41.stan

data{
  int N;
  array[N] vector[2] y;
}

parameters{
  vector[2] m;
  vector<lower=0>[2] s;
  real<lower=-1,upper=1> r;
}

transformed parameters{
  cov_matrix[2] cv;
  cv[1,1]=s[1]^2;
  cv[2,2]=s[2]^2;
  cv[1,2]=r*s[1]*s[2];
  cv[2,1]=r*s[1]*s[2];
}

model{
  y~multi_normal(m,cv);
}

ex4-42.stan

data{
  int N;
  matrix[2,2] dp;
}

parameters{
  vector<lower=0>[2] s;
  real<lower=-1,upper=1> r;
}

transformed parameters{
  cov_matrix[2] cv;
  cv[1,1]=s[1]^2;
  cv[2,2]=s[2]^2;
  cv[1,2]=r*s[1]*s[2];
  cv[2,1]=r*s[1]*s[2];
}

model{
  dp~wishart(N-1,cv);
}

ex4-43.stan

data{
  int N;
  int M;
  matrix[M,M] dp;
}

parameters{
  vector<lower=0>[M] s;
  corr_matrix[M] r;
}

transformed parameters{
  cov_matrix[M] cv;
  cv=quad_form_diag(r,s);
}

model{
  dp~wishart(N-1,cv);
}

ex4-44.stan

data{
  int N;
  int K;
  array[N] vector[K] y;
}
parameters{
  vector[K] m;
  cov_matrix[K] cov;
}
model{
  y~multi_normal(m,cov);
}
library(MASS)
n=20
y=mvrnorm(n,c(0,0),matrix(c(1,0.7,0.7,1),2),2)
cov(y)
##       [,1]  [,2]
## [1,] 0.867 0.538
## [2,] 0.538 0.784
cor(y)
##       [,1]  [,2]
## [1,] 1.000 0.653
## [2,] 0.653 1.000
plot(y)

#estimate covariance
mdl=cmdstan_model('./ex4-41.stan')

data=list(N=n,y=y)

mle=mdl$optimize(data=data)
## Initial log joint probability = -35.5849 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       17      -9.54568   2.07811e-05   0.000529216      0.7427      0.7427       18    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##   lp__       -9.55
##   m[1]        0.07
##   m[2]        0.25
##   s[1]        0.91
##   s[2]        0.86
##   r           0.65
##   cv[1,1]     0.82
##   cv[2,1]     0.51
##   cv[1,2]     0.51
##   cv[2,2]     0.74
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,ch=F)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##   lp__    -13.71 -13.37 1.72 1.48 -16.90 -11.59 1.01      825     1159
##   m[1]      0.07   0.07 0.24 0.23  -0.32   0.45 1.00     1295     1239
##   m[2]      0.24   0.25 0.24 0.22  -0.14   0.62 1.00     1355     1330
##   s[1]      1.01   0.98 0.18 0.16   0.76   1.34 1.00     1284     1300
##   s[2]      0.96   0.93 0.18 0.16   0.71   1.30 1.00     1183     1251
##   r         0.61   0.63 0.15 0.14   0.32   0.81 1.02      443      679
##   cv[1,1]   1.05   0.97 0.40 0.31   0.58   1.79 1.00     1284     1300
##   cv[2,1]   0.62   0.55 0.31 0.26   0.24   1.21 1.01      550      644
##   cv[1,2]   0.62   0.55 0.31 0.26   0.24   1.21 1.01      550      644
##   cv[2,2]   0.96   0.87 0.39 0.29   0.51   1.69 1.00     1183     1251

#estimate correlation,covariance
#deviation product sum matirix folows Wishart dist.
#deviation product sum = covariance * n-1 or n

mdl=cmdstan_model('./ex4-42.stan')

data=list(N=n,dp=cov(y)*(n-1))

mle=mdl$optimize(data=data)
## Initial log joint probability = -324.231 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       14       -10.043   9.20225e-05   0.000306709           1           1       18    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##   lp__      -10.04
##   s[1]        0.93
##   s[2]        0.89
##   r           0.65
##   cv[1,1]     0.87
##   cv[2,1]     0.54
##   cv[1,2]     0.54
##   cv[2,2]     0.78
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,ch=F)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##   lp__    -13.04 -12.72 1.32 1.05 -15.60 -11.63 1.01      830      687
##   s[1]      1.01   0.99 0.19 0.17   0.76   1.34 1.00     1395     1240
##   s[2]      0.95   0.93 0.17 0.16   0.72   1.26 1.00     1381     1160
##   r         0.60   0.62 0.16 0.14   0.30   0.80 1.01      455      506
##   cv[1,1]   1.06   0.98 0.44 0.33   0.58   1.80 1.00     1395     1240
##   cv[2,1]   0.60   0.55 0.31 0.25   0.23   1.14 1.01      564      603
##   cv[1,2]   0.60   0.55 0.31 0.25   0.23   1.14 1.01      564      603
##   cv[2,2]   0.94   0.87 0.36 0.30   0.52   1.59 1.00     1381     1160

#estimate correlation,covariance
mdl=cmdstan_model('./ex4-43.stan')

m=2
data=list(N=n,M=m,dp=cov(y)*(n-1))

mle=mdl$optimize(data=data)
## Initial log joint probability = -1230.72 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       16       -10.043   5.06124e-05   1.12298e-05           1           1       21    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##   lp__      -10.04
##   s[1]        0.93
##   s[2]        0.89
##   r[1,1]      1.00
##   r[2,1]      0.65
##   r[1,2]      0.65
##   r[2,2]      1.00
##   cv[1,1]     0.87
##   cv[2,1]     0.54
##   cv[1,2]     0.54
##   cv[2,2]     0.78
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,ch=F)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##   lp__    -12.30 -11.97 1.28 1.04 -14.79 -10.92 1.00      812     1252
##   s[1]      1.01   0.98 0.18 0.17   0.76   1.34 1.00      953     1146
##   s[2]      0.96   0.94 0.17 0.16   0.73   1.27 1.00     1121     1241
##   r[1,1]    1.00   1.00 0.00 0.00   1.00   1.00   NA       NA       NA
##   r[2,1]    0.61   0.63 0.15 0.14   0.33   0.81 1.00      735     1032
##   r[1,2]    0.61   0.63 0.15 0.14   0.33   0.81 1.00      735     1032
##   r[2,2]    1.00   1.00 0.00 0.00   1.00   1.00   NA       NA       NA
##   cv[1,1]   1.05   0.96 0.41 0.33   0.58   1.80 1.00      953     1146
##   cv[2,1]   0.62   0.56 0.30 0.25   0.25   1.20 1.00      684      747
##   cv[1,2]   0.62   0.56 0.30 0.25   0.25   1.20 1.00      684      747
##   cv[2,2]   0.95   0.88 0.37 0.30   0.53   1.63 1.00     1121     1241

mdl=cmdstan_model('./ex4-44.stan')

k=2
data=list(N=n,K=k,y=y)

mle=mdl$optimize(data=data)
## Initial log joint probability = -1110.98 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       17      -9.54568   1.35409e-05   0.000670983      0.5952      0.5952       24    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##  lp__        -9.55
##  m[1]         0.07
##  m[2]         0.25
##  cov[1,1]     0.82
##  cov[2,1]     0.51
##  cov[1,2]     0.51
##  cov[2,2]     0.74
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,ch=F)
##  variable   mean median   sd  mad     q5   q95 rhat ess_bulk ess_tail
##  lp__     -12.04 -11.70 1.77 1.67 -15.39 -9.82 1.01      827     1236
##  m[1]       0.08   0.07 0.26 0.24  -0.35  0.51 1.00     1062     1160
##  m[2]       0.25   0.25 0.24 0.22  -0.14  0.64 1.00     1095     1033
##  cov[1,1]   1.27   1.14 0.55 0.42   0.65  2.35 1.00     1330     1338
##  cov[2,1]   0.78   0.69 0.44 0.35   0.26  1.61 1.00      926      891
##  cov[1,2]   0.78   0.69 0.44 0.35   0.26  1.61 1.00      926      891
##  cov[2,2]   1.14   1.03 0.50 0.39   0.56  2.07 1.00     1209     1127